1. Design of Metamaterials using Topology Optimization

Property-based topology optimization designs a microstructure whose homogenized stiffness matches a target tensor. Use it to synthesize metamaterial unit cells with SIMP, filters and sensitivity weighting, then to verify the result with a full-scale finite-element model.

Reading time

9 min

Length

1,577 words

Updated

Sep 15, 2026

Total views

--

#mechanics/metamaterial#mechanics/topology-optimization

1. Topology Optimization using the target mechanical property

There are two main reference for the property-based topology optimization, including ^[1] and ^[2]

We already know the homogenization issue and the energy-based homogenization from ^[3] and ^[4]. So we start from the Section4 of ^[1:1], which is the optimization model.

(1) Optimization Model

We use the SIMP approach for the design of metamaterials, the density is defined as :

$$\begin{array}{}\text{(1.1.1)}& E_e({\rho }_e)=E_{min}+{\rho }^p(E_0-E_{min})\end{array}$$

The Formulation of optimization should be :

$$\begin{array}{}\text{(1.1.2)}& \begin{array}{rl}& \underset{\rho }{min}:c(E_{ijkl}^H)(\rho )k,l=1,\dots ,d\\ & \text{s.t.}:K{U}^{A(kl)}=F^{(kl)}\\ & \sum _{e=1}^N{v}_e\frac{{\rho }_e}{|Y|}\le \theta \\ & 0\le {\rho }_e\le 1\end{array}\end{array}$$

Here, $c$ is an extractor function for the $E_{ijkl}^H$. Where for 2D case, $i,j,k,l=1,2$.

(2) FE Analysis Process

Firstly, to check how periodic boundary conditions are applied and solved, for example MATLAB code of it, see ^[5] for details.

For the FE analysis, we can use follows :

$$\begin{array}{}\text{(1.2.1)}& \left[\begin{array}{cccc}{K}_{11}& K_{12}& \dots & K_{14}\\ \dots \\ \dots \\ & & & K_{44}\end{array}\right]\left[\begin{array}{c}{\bar{U}}_1\\ U_2\\ U_3\\ U_4\end{array}\right]=\left[\begin{array}{c}{F}_1\\ F_2\\ F_3\\ F_4\end{array}\right]\end{array}$$

We note $F_1$ is unknown, $F_2=0$, $F_3+F_4=0$, then $U_4=U_3+\bar{W}$ and $K$ is symmetric

Reduction to the $2\times 2$ system (Matrix reduction technique)

The prescribed displacement and periodic relation leave $U_2$ and $U_3$ as the independent unknowns. We should define the unknown vector $q=[U_2,U_3{]}^T$, and split form of $U$

$$\begin{array}{}\text{(1.2.2)}& \mathbf{q}=\left[\begin{array}{c}{U}_2\\ U_3\end{array}\right],{\mathbf{U}}_0=\left[\begin{array}{c}{\bar{U}}_1\\ 0\\ 0\\ \bar{W}\end{array}\right],\mathbf{T}=\left[\begin{array}{cc}0& 0\\ 1& 0\\ 0& 1\\ 0& 1\end{array}\right].\end{array}$$

Therefore, the complete displacement vector is :

$$\begin{array}{}\text{(1.2.3)}& \mathbf{U}={\mathbf{U}}_0+\mathbf{T}\mathbf{q}.\end{array}$$

Premultiplication by ${\mathbf{T}}^T$ removes the unknown reaction $F_1$ and combines the two periodic-boundary equations into $F_3+F_4=0$ :

$$\begin{array}{}\text{(1.2.4)}& {\mathbf{T}}^T\mathbf{K}\mathbf{T}\mathbf{q}=-{\mathbf{T}}^T\mathbf{K}{\mathbf{U}}_0.\end{array}$$

Using the symmetry of $\mathbf{K}$ gives the explicit reduced system

$$\begin{array}{}\text{(1.2.5)}& \left[\begin{array}{cc}{K}_{22}& K_{23}+K_{24}\\ K_{23}+K_{24}& K_{33}+2K_{34}+K_{44}\end{array}\right]\left[\begin{array}{c}{U}_2\\ U_3\end{array}\right]=-\left[\begin{array}{c}{K}_{12}{\bar{U}}_1+K_{24}\bar{W}\\ (K_{13}+K_{14}){\bar{U}}_1+(K_{34}+K_{44})\bar{W}\end{array}\right].\end{array}$$

After solving this system, recover the slave-boundary displacement and the reaction by

$$\begin{array}{}\text{(1.2.6)}& U_4=U_3+\bar{W},F_1=(\mathbf{K}\mathbf{U}{)}_1.\end{array}$$

Matrix Reduction
For a general matrix with some known variables, we can use following steps to reduce it to a lower-rank problem :

$$Ax=b$$

we don't want $b$ to appear in the final form.

1. If the $x$ is partially known, express the $x$ as :

$$x=x_0+Tq$$

where $q$ is unknown vectors

2. Then we need another matrix $P$

$$P^{T}b=0\to P^{T}K(x_0+Tq)=0$$

the new form is reached.

This transformation is also directly applicable in symbolic software such as SymPy, Mathematica, or the MATLAB Symbolic Math Toolbox. In an FE solver, the same relation is generally imposed with a multi-point constraint (MPC).

%%%%%%%%%%%%%%%%%% FE-ANALYSIS  %%%%%%%%%%%%%%
sK = reshape(KE(:)*(Emin+xPhys(:)'.^penal*(E0-Emin)),64*nelx*nely,1);
K = sparse(iK,jK,sK); K = (K+K')/2;
Kr = [K(d2,d2), K(d2,d3)+K(d2,d4); K(d3,d2)+K(d4,d2), K(d3,d3)+K(d4,d3)+K(d3,d4)+K(d4,d4)];
U(d1,:) = ufixed;
U([d2,d3],:) = Kr\(-[K(d2,d1); K(d3,d1)+K(d4,d1)]*ufixed-[K(d2,d4); K(d3,d4)+K(d4,d4)]*wfixed);
U(d4,:) = U(d3,:)+wfixed;

(3) Derivatives of Homogenized Stiffness matrix

We compute the derivation against every element density ${\rho }_e$ as :

$$\begin{array}{}\text{(1.3.1)}& \frac{\partial E_{ijkl}^H}{\partial {\rho }_e}=\frac{1}{|Y|}p{\rho }_e^{p-1}\end{array}$$

From ^[2:1], we know that we have (1.1.1)^[3:1], then :

$$\begin{array}{}\text{(1.3.2)}& E_{ijkl}^H=\frac{1}{|Y|}\sum _{e=1}^{N_E}(u_e^0-u_e^{\ast (ij)})k_e(u_e^0-u_e^{\ast (kl)})\end{array}$$

Then we have $k_e=\rho k_e^0$, then the derivation is :

$$\begin{array}{}\text{(1.3.3)}& \boxed{{\displaystyle \frac{\partial E_{ijkl}^H}{\partial {\rho }_e}=\frac{1}{Y}p{\rho }^{p-1}(u_e^0-u_e^{\ast (ij)})k_e(u_e^0-u_e^{\ast (kl)})}}\end{array}$$

To compute above expression in dolfinx, it's an anisotropic computation, the stress of it should be computed by $C_{66}$. We note that in the computation of homogenization, we have already solved $u_e^{ij}$, so we can use mat.mu and mat.lmbda to compute the stress :

$$\begin{array}{}\text{(1.3.4)}& ϵ={ϵ}^{\sigma }+{ϵ}^{\text{th}}+{ϵ}^{\text{v}}\end{array}$$

(4) Application cases for the design of meta-materials

For example, we want to maximize the bulk modulus^[6] in (1.7)^[6:1], we will minimize the negative of it :

$$\begin{array}{}\text{(1.4.1)}& c(E_{ijkl})=-4K_{\epsilon }^H(\rho )\stackrel{2D}{⟶}-(E_{1111}+E_{1122}+E_{2211}+E_{2222})\end{array}$$

Or the shear modulus in (2.4)^[6:2], it is :

$$\begin{array}{}\text{(1.4.2)}& c(E_{ijkl})=-E_{ijij}\stackrel{2D}{⟶}-E_{1212}\end{array}$$

the $i,j$ ==depends on the direction it wants to minimize==.

We also note besides (1.4.1~2), the article^[2:2] gives a more general form, reads :

$$\begin{array}{}\text{(1.4.3)}& \text{Maximize :}J(u,\mathrm{\Phi })=\sum _{i,j,k,l=1}^d{\eta }_{ijkl}{E}_{ijkl}^H(\rho )\end{array}$$

where ${\eta }_{ijkl}$ denotes the coefficient of each term.

If we want to reach a material with negative-poission's ratio, we want to minimize the Poisson's ratio $\nu$ , Using (2.3)^[7] and Fig. 2.2^[7:1], we have :

$$\begin{array}{}\text{(1.4.4)}& \cancel{\nu =\frac{E_{1122}}{E_{1111}+E_{1122}}}\end{array}$$

The above expression is not the general Poisson ratio of a homogenized 2D cell. It recovers the original 3D isotropic Poisson ratio only when (C) is an isotropic plane-strain stiffness tensor.

Since the optimized microstructure is generally anisotropic, so that conversion no longer applies.

The definition under uniaxial loading in direction 1 is :

$$\begin{array}{}\text{(1.4.5)}& {\nu }_{12}-\frac{{\epsilon }_{22}}{{\epsilon }_{11}}\text{when}{\sigma }_{22}={\sigma }_{12}=0\end{array}$$

We let the homogenized compliance as :

$$\begin{array}{}\text{(1.4.6)}& S^H=(C^H{)}^{-1}\end{array}$$

Then the general expression of Poisson's ratio are :

$$\begin{array}{}\text{(1.4.7)}& {\nu }_{12}=-\frac{S_{2211}^H}{S_{1111}^H}\end{array}$$

and :

$$\begin{array}{}\text{(1.4.8)}& {\nu }_{21}=-\frac{S_{1122}^H}{S_{2222}^H}\end{array}$$

To get the poisson's ratio, we can use :

C_sym = 0.5 * (C_H + C_H.T)
S_H = np.linalg.inv(C_sym)

In the original essay, the target function is set as:

$$\begin{array}{}\text{(1.4.9)}& c=E_{1122}-{\beta }^l(E_{1111}+E_{2222})\end{array}$$

where $l$ is the loop number, which finally transfer the objective to elastic tensor $E_{1122}$. And also, to prevent the OC issue, the essay modified the OC bisection rule as :

xnew = max(0,max(x-move,min(1, min(x+move,x.*(-dc./dv/lmid)))));

Why setting this target function reasonable ? Under uniaxial stress in direction 1, we have :

$$S=\frac{C^{\ast }}{|C|}\to v_{12}\to \frac{C_{1122}}{C_{2222}}$$

To maximize $C_{1122}$ while minimize $C_{2222}$ and $C_{1111}$, the target $c$ is reasonable.

2. Implementation by FEniCSx

(1) Indexing rule of Sphinx

1. In the dolfinx, every component has a ufl_shape attribute, so they have multiple components, We can access the components directly by indexing. For example, 2d strain/stress ${\epsilon }_{ij}$ is (2x2) matrix and 3d strain is (3x3) matrix.

u_sol[0]   # x displacement  
u_sol[1]   # y displacement 

# they returns a Indexed object 
# For example, to compute the magnitude in another space, we can use interpolate + expression : 
V_mag = fem.functionspace(msh, ("CG", 1))
u_magnitude = fem.Function(V_mag, name="displacement_magnitude")
u_magnitude.interpolate(
    fem.Expression(
        ufl.sqrt(u_sol[0] ** 2 + u_sol[1] ** 2),  # also ufl.inner(u_sol, u_sol)
        V_mag.element.interpolation_points,
    )
)
u_magnitude.x.scatter_forward()
plot_function_pv(u_magnitude, title="displacement magnitude")

(2) Stress Computation for anisotropic materials

So, having the C66 matrix $C^H$, we also need the voigt strain form, which is :

$$\begin{array}{}\text{(2.2.1)}& {\left[\begin{array}{ccc}{\epsilon }_{11}& {\epsilon }_{22}& {\epsilon }_{12}\end{array}\right]}^T(2D){\left[\begin{array}{cccccc}{\epsilon }_{11}& {\epsilon }_{22}& {\epsilon }_{22}& {\epsilon }_{12}& {\epsilon }_{13}& {\epsilon }_{23}\end{array}\right]}^T(3D)\end{array}$$

In implementation, we firstly, need two transfer functions to transfer the tensor notation to voigt notation, as given in ^[8], then, we can implement the stress computation as :

def stress_e_c66 (C66, eps):
	eps_voigt = tensor_to_voigt(eps, engineering_shear=True)
	C_const = ufl.as_matrix(C) 
	sigma_voigt = ufl.dot(
        C_ufl,
	    eps_voigt,
    )
	return voigt_to_tensor(sigma_voigt, dim)

(3) Computation of Derivatives of Stiffness Matrix

We know from (1.3.3) that, we can compute the derivatives by :

$$\begin{array}{}\text{(2.3.1)}& \frac{\partial E_{ijkl}^H}{\partial {\rho }_e}=p{\rho }^{p-1}({\epsilon }_e^0-{\epsilon }_e^{ij})C({\epsilon }_e^0-{\epsilon }_e^{kl})\end{array}$$

For sensitivity of stiffness matrix, we have following computations :

def effective_stiffness_sensitivity(
    homo: CellHomogenization,
    rho: np.ndarray,
) -> tuple[np.ndarray, np.ndarray]:
    r"""Assemble the elementwise derivatives of the homogenized stiffness.

    .. math::
        \frac{\partial C^H_{IJ}}{\partial \rho_e}
        = \frac{1}{\lvert Y \rvert} \int_{Y_e} \eta^{(I)} : \mathbb{C}_{,\rho} : \eta^{(J)} \, dY.

    where :math:`I = ij,  J = kl`

    The EBHM correctors use
    :math:`\eta^{(I)} = E^{(I)} - \varepsilon(\chi^{(I)})`.
    """
    homo.update_density(rho)
    mat, msh = homo.mat, homo.msh
    V0 = fem.functionspace(msh, ("DG", 0))
    C_H, correctors = homo.solve_ebhm(return_u_solu=True)
    dim = homo.dim
    comps = len(unit_strain_tensors(dim=dim))
    if len(correctors) != comps:
        raise ValueError("Number of correctors must match the Voigt components.")

    unit_strains = unit_strain_tensors(dim=dim)
    eta = [
        fem.Constant(msh, unit_strains[I]) - strain(correctors[I])
        for I in range(comps)
    ]
    unit_vec = ufl.TestFunction(V0)
    dC_H_drho = np.zeros((comps, comps, len(rho)))

    dscale = (
        homo.penalty
        * (homo.rho_max - homo.rho_min)
        * rho ** (homo.penalty - 1.0)
    )
    for I, eta_I in enumerate(eta):
        for J, eta_J in enumerate(eta):
            Q = ufl.inner(stress_e(mat.mu, mat.lmbda, eta_I), eta_J)
            Q_vec = fem.assemble_vector(
                fem.form(Q * unit_vec * ufl.dx)
            )
            dC_H_drho[I, J, :] = dscale * Q_vec.array / homo.vol
    return C_H, dC_H_drho

We note here we use the bulk modulus in (1.4.1) as the objective, After getting the derivative, we can plot and check the components, the [0, 0] component of the gradient $\frac{\partial C}{\partial \rho }$ is plotted as follows :

while the [0, 1] component is like :

(4) Essential problems

1) Filter Configurations

Firstly we have, for (2.3.1) :

$$\begin{array}{}\text{(2.4.1)}& \frac{\partial C_{IJ}^H}{\partial {\rho }_e}=\frac{1}{|Y|}{\int }_{Y_e}p({\rho }_{max}-{\rho }_{min}){\rho }_e^{p-1}{\eta }^{(I)}:{\mathbb{C}}_0:{\eta }^{(J)}dY\end{array}$$

For the density filtering, we only use the density filtering only to get the final physical density, if we apply the filter in the middle process, the filter may hinder the density transformation process, or even deteriorate the results.

If we use the physical (filtered) density to compute the $\frac{dJ}{d\rho }$, the result is $\frac{\partial J}{\partial {\rho }_p}$. so we have $\delta {\rho }_p=H\delta {\rho }_d$ where $p$ means physical and $d$ means design :

$$\begin{array}{}\text{(2.4.2)}& \delta J={\left(\frac{\partial J}{\partial {\rho }_p}\right)}^T\delta {\rho }_p=\underset{\text{filtered density }}{\underbrace{{\left(\frac{\partial J}{\partial {\rho }_p}\right)}^{T}H}}\delta {\rho }_d\end{array}$$

We have following result :

$$\begin{array}{}\text{(2.4.3)}& \boxed{{\displaystyle \frac{\partial J}{\partial \rho }{|}_{\text{filtered}}=H^T\frac{\partial J}{\partial {\rho }_p}}}\end{array}$$

In the code example, we can implement this filtering process by following :

density_filter = DistanceBasedFilter(  
    msh,  
    rmin=filter_radius * min_mesh_size(msh),  
    p=1.0,  
    filter_type="element",  
)
C_H, dC_H_drho = effective_stiffness_sensitivity(
	homo,
	rho_physical,
)
# C^H_1111 + C^H_1122 + C^H_2211 + C^H_2222
bulk_stiffness = float(np.sum(C_H[:2, :2]))  # float(C_H[0, 1])  
dJ_drho = -np.sum(dC_H_drho[:2, :2, :], axis=(0, 1))  # - dC_H_drho[0, 1, :]    
dJ_drho_filtered = density_filter.apply_adjoint(dJ_drho)
vols_filtered = density_filter.apply_adjoint(vols)
Lmbda = optimization_criteria(
	rho=rho_design,
	dc_drho=dJ_drho_filtered,
	elem_vols=vols_filtered,
	vol_frac=volume_fraction,
	zeta=move_limit,
	precision=1e-6,
	lmbda_range=(1e-30, 1e30),
	rho_min=rho_min,
	rho_max=rho_max,
)

In the essay that we references, the article provides 2 types of filters. In our cellgen package, the DistanceBasedFilter applying at both sensitivity and volume is same to filter2. According to the article, there are 2 types of filters, given in the following code :

%% FILTERING/MODIFICATION OF SENSITIVITIES
if ft == 1
dc(:) = H*(x(:).*dc(:))./Hs./max(1e-3,x(:));
elseif ft == 2
dc(:) = H*(dc(:)./Hs);
dv(:) = H*(dv(:)./Hs);
end

$$\frac{\partial J}{\partial \rho }{|}_{\text{filtered}}=\frac{1}{\rho }{H}^T\left(\rho \frac{\partial J}{\partial \rho }\right)$$

2) The crucial Initial Condition Problem

In the debugging process, one of the critical problem is the initial condition. This can simply be fixed by modifying the line as :

if sqrt((i-nelx/2-0.5)^2+(j-nely/2-0.5)^2) < min(nelx,nely)/3 
%%%%%%%%%%%% set to follows %%%%%%%%%%%%
if sqrt((i-nelx/2-0.5)^2+(j-nely/2-0.5)^2) < 0.15 * nelx

3) MMA problem for inner density increase

The unconservative MMA updates, while topX.m uses density-filtered OC bisection. MMA continued making large moves after the objective peaked, transferring density into the center. So, the OC optimizer is generally stabler than MMA optimizer.

We note that lower penalty is also beneficial. penalty = 3.5 -> 4 is a better choice than the penalty = 5.

One of the problem is, the filter radius is also crucial. A good choice is $1.3\sim 1.5$. When the filter radius is set to a larger value like 2, the result can be asymmetric. The following structure is reached by setting r=1.5 .

The following gives a result on $r=2$ :

(5) Poisson's ratio

By the optimization of Poisson's ratio, we can still reach some materials with negative Poisson's ratio.

The generation configs is center hole radius 1/6, with following configs :

volume_fraction = 0.3
rho_min = 1e-9

config = UnitCellOptimProblemConfig(
	vol_frac=volume_fraction,
	max_iterations=150,
	result_dir=str(Path(__file__).parent),
	result_filename="microcell_top_design_density.xdmf",
	density_filename="microcell_top_design_best_density.npy",
	penalty=3.5,
	filter_radius=1.2,
	rho_min=rho_min,
	rho_max=1.0,
	move_limit=0.1,
	backend="mma",  # Use "oc" to select the optimality-criteria update.
	mma_auxiliary_penalty=1e5,
	mma_max_subproblem_iter=200,
)

Under these configurations, the poission's ratio is -0.5, resulting microstructure is given in Fig. 2.7 :

We can validate the auxetic properties by ANSYS using a lattice constructed by this cell :

Another typical configuration is:

def center_hole(x):
    return (x[0, :] - 0.5) ** 2 + (x[1, :] - 0.5) ** 2 < ((0.2) ** 2)

volume_fraction = 0.4
rho_min = 1e-9

config = UnitCellOptimProblemConfig(
	vol_frac=volume_fraction,
	max_iterations=150,
	result_dir=str(Path(__file__).parent),
	result_filename="microcell_top_design_density.xdmf",
	density_filename="microcell_top_design_best_density.npy",
	penalty=3.5,
	filter_p=0.7,
	filter_radius=1.0,
	rho_min=rho_min,
	rho_max=1.0,
	move_limit=0.1,
	backend="mma",  # Use "oc" to select the optimality-criteria update.
	mma_auxiliary_penalty=1e5,
	mma_max_subproblem_iter=200,
)

The optimization result under 50x50 resolutions are given in Fig. 2.9

We note fhat we may need to increase the resolution into 200x200 for a detailed structure, We can get the structure with more detailed structure layout, given by following picture :

---

1. Xia, Liang, and Piotr Breitkopf. “Design of Materials Using Topology Optimization and Energy-Based Homogenization Approach in Matlab.” Structural and Multidisciplinary Optimization 52, no. 6 (2015): 1229–41. https://doi.org/10.1007/s00158-015-1294-0. ↩︎ ↩︎

2. Gao, Jie, Hao Li, Zhen Luo, Liang Gao, and Al Et. “Topology Optimization of Micro-Structured Materials Featured with the Specific Mechanical Properties.” International Journal of Computational Methods, ahead of print, 2019. https://doi.org/10.1142/S021987621850144X. ↩︎ ↩︎ ↩︎

3. Second Form of Equivalent Stiffness ↩︎ ↩︎

4. EBHM ↩︎

5. PBC ↩︎

6. 3. Definition of Common Mechanical Modulus ↩︎ ↩︎ ↩︎

7. Hook's Law ↩︎ ↩︎

8. 1. Symmetric Stiffness Tensor and Voigt notation ↩︎